Control method, controller, recording medium recording control program, numerical calculation method, numerial calculator and recording medium recording numerical calculation program

ABSTRACT

A·U+B(U)=f, wherein A is a linear differential operator, B is a nonlinear differential operator, and f is an inhomogeneous term (source term) in a nonlinear partial differential equation to be satisfied by a physical quantity U, is solved by successive approximation. In calculation, (f−A·U m −B(U m )) is given as a nonlinear residual rr of an approximate solution U m , wherein m is the number of repeating times, and the approximate solution U m  is repeatedly corrected so as to reduce a norm of a nonlinear residual r m+1  employed in a subsequent step.

BACKGROUND OF THE INVENTION

[0001] The present invention relates to a technique to control anonlinear operation for a physical quantity such as a pressure or atemperature, and more particularly, it relates to a control techniqueusing successive approximation algorithm. Also, it relates to anumerical calculation technique for a physical quantity.

[0002] In the technique to control a nonlinear operation for a physicalquantity such as a pressure or a temperature, successive approximationalgorithm has been conventionally used. In the successive approximationalgorithm, successive approximation is repeated on the basis of avariety of optimization theories so as to obtain an optimum solutionthrough numerical calculation.

[0003] In the conventional successive approximation algorithm, however,it is premised that a nonlinear operation can be expressed byaccumulating small linear operations obtained by linearizing a nonlinearequation. Therefore, in the case where a phenomenon with strongnonlinearity occurs, the calculation should be repeated a massive numberof times for obtaining solutions of respective parameters used fordetermining the control operation. Accordingly, rapid and highly precisecontrol is disadvantageously difficult due to the time spent on thecalculation, or the calculation is diverged while it is repeated andhence a solution cannot be obtained, which disables the controloperation.

SUMMARY OF THE INVENTION

[0004] An object of the invention is providing a control technique and anumerical calculation technique using successive approximation algorithmfor controlling a nonlinear operation more rapidly and more preciselythan in the conventional technique.

[0005] Specifically, according to the invention, as a technique tocontrol or numerically calculate a physical quantity U, A·U+B(U)=f,wherein A is a linear differential operator, B is a nonlineardifferential operator, and f is an inhomogeneous term (source term) in anonlinear partial differential equation to be satisfied by the physicalquantity U, is solved by using successive approximation. In thecalculation, (f−A·U^(m)−B(U^(m))) is given as a nonlinear residual r^(m)of an approximate solution U^(m), wherein m is the number of repeatingtimes, and the approximate solution U^(m) is repeatedly corrected so asto reduce a nonlinear residual r^(m+1) employed in a subsequent step.

[0006] In the control or numerical calculation technique of thisinvention, the nonlinear term B(U) is expressed as follows:${B(U)} = {\sum\limits_{n}{{\overset{n}{N}(U)}{\overset{n}{D} \cdot U}}}$

${\overset{n}{N}(U)}\quad \cdots$

[0007] . . . nonlinear coefficient, ${\overset{n}{D}(U)}\quad \cdots$

[0008] . . . differential operator N(U^(m)+φ) is expanded with respectto the approximate solution U^(m) and a perturbation quantity φ by theTaylor expansion or the mean value theorem, to give the following withterms higher than φ² ignored: $\begin{matrix}{{{\left. {{B\left( U^{m} \right)} + \varphi} \right) \cong {{B\left( U^{m} \right)} + {L \cdot \varphi} + {\varphi {\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot \varphi}}}}}},\quad \ldots}{{{L \cdot \varphi} \equiv {\sum\limits_{n}{\left\lbrack {{{\overset{n}{N}\left( U^{m} \right)}\overset{n}{D}} + {{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot U^{m}}}} \right\rbrack \varphi}}},{{\overset{n}{N^{\prime}}\left( U^{m} \right)} \equiv \left\lbrack \frac{\partial{\overset{n}{N}(U)}}{\partial U} \right\rbrack_{U = U^{m}}}}} & (4)\end{matrix}$

[0009] and the following first process and second process are preferablyrepeatedly executed until the approximate solution U^(m) is converged.The first process includes the steps of setting U⁰ as an initial valueof the physical quantity U; setting 0 as an initial value of the numberm of repeating times and (f−A·U⁰−B(U⁰)) as an initial value r⁰ of anonlinear residual r; and obtaining a predicted approximate value φ^(m)of the following equation through repeated calculations whileincrementing the number m of repeating times: $\begin{matrix}{{\left\lbrack {A + L + {\varphi {\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}\overset{n}{D}}}}} \right\rbrack \cdot \varphi} = r^{m}} & (5)\end{matrix}$

[0010] The second process includes the steps of obtaining a correctedapproximate value φ^(m) for minimizing a norm of a nonlinear residualr^(m+1) on the basis of a predicted approximate values φ^(m) andcorrected approximate values φ^(m−1), . . . , and φ^(m−Lmax+1), whereinLmax is an integer of 2 or more; giving (U^(m)+φ^(m)) as an approximatesolution U^(m+1); and giving (f−A·U^(m+1)−B(U^(m+1))) as the nonlinearresidual r^(m+1).

[0011] Alternatively, in the control or numerical calculation techniqueof this invention, the nonlinear term B(U) is expressed as follows:${B(U)} = {\sum\limits_{n}{{\overset{n}{N}(U)}{\overset{n}{D} \cdot U}}}$

${\overset{n}{N}(U)}\quad \cdots$

[0012] . . . nonlinear coefficient, ${\overset{n}{D}(U)}\quad \cdots$

[0013] . . . differential operator N(U^(m)+φ)) is expanded with respectto the approximate solution U^(m) and a perturbation quantity φ by theTaylor expansion or the mean value theorem, to give the following withterms higher than φ ignored:

B(U ^(m)+φ)≅B(U ^(m))+L·φ,  (4′)${{L \cdot \varphi} \equiv {\sum\limits_{n}{\left\lbrack {{{\overset{n}{N}\left( U^{m} \right)}\overset{n}{D}} + {{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot U^{m}}}} \right\rbrack \varphi}}},{{\overset{n}{N^{\prime}}\left( U^{m} \right)} \equiv \left\lbrack \frac{\partial{\overset{n}{N}(U)}}{\partial U} \right\rbrack_{U = U^{m}}}$

[0014] and the following first process and second process are preferablyrepeatedly executed until the approximate solution U^(m) is converged.The first process includes the steps of setting U⁰ as an initial valueof the physical quantity U; setting 0 as an initial value of the numberm of repeating times and (f−A·U⁰−B(U⁰)) as an initial value r⁰ of anonlinear residual r; obtaining a predicted approximate value φ^(m) of[A+L]φ=r^(m) through repeated calculations while incrementing the numberm of repeating times. The second process includes the steps of obtaininga corrected approximate value φ^(m) for minimizing a norm of a nonlinearresidual r^(m+1) on the basis of a predicted approximate value φ^(m) andcorrected approximate values φ^(m−1), . . . and φ^(m−Lmax+1), whereinLmax is an integer of 2 or more; giving (U^(m)+φ^(m)) as an approximatesolution U^(m+1); and giving (f−A·U^(m+1)−B(U^(m+1))) as the nonlinearresidual r^(m+1).

BRIEF DESCRIPTION OF THE DRAWINGS

[0015]FIG. 1 is a block diagram for showing the architecture of acontroller according to an embodiment of the invention;

[0016]FIG. 2 is a flowchart for showing procedures in a control methodaccording to Embodiment 1 of the invention;

[0017]FIG. 3 is a flowchart for showing procedures in a control methodaccording to Embodiment 2 of the invention;

[0018]FIGS. 4A and 4B are diagrams for showing an example of gridgeneration to which algorithm of this invention is applied;

[0019]FIGS. 5A and 5B are graphs for comparing the convergence rate ofthe example shown in FIGS. 4A and 4B;

[0020]FIGS. 6A and 6B are diagrams for showing another example of thegrid generation to which the algorithm of this invention is applied; and

[0021]FIGS. 7A and 7B are graphs for comparing the convergence rate ofthe example shown in FIGS. 6A and 6B.

DETAILED DESCRIPTION OF THE INVENTION

[0022] Preferred embodiments of the invention will now be described withreference to the accompanying drawings.

[0023]FIG. 1 is a block diagram for showing the architecture of acontroller according to an embodiment of the invention, and thisarchitecture is shown as an example of application of a controltechnique according to this invention to vehicle automatic drivecontrol. It is assumed in FIG. 1 that a vehicle is equipped with acircuit for allowing the vehicle to drive to a destination along anorbit set by a destination/route setting unit 11 within an allowableerror range. During the drive, drive information such asposition/azimuth information and a driving speed are detected by aposition information detector 12. Also, information such as a set valueand a driving state are displayed in a real-time manner on a display 18.

[0024] A shift φ between the position of the vehicle and the orbit isalways monitored by information management unit 13. If the shift φexceeds an allowable value δ_(max), by using the drive informationobtained at this point as starting values, optimization parametercontrol of an azimuth angle, angular acceleration, a speed, accelerationand the like, namely, optimization means for returning the vehicle tothe orbit, is simulated by using algorithm of this invention.

[0025] As a characteristic of this invention, a nonlinear residual isdirectly used in repeated calculations of the successive approximationalgorithm.

[0026] Embodiment 1

[0027]FIG. 2 is a flowchart for showing procedures in a control methodaccording to Embodiment 1 of the invention.

[0028] <Formulation of Nonlinear Residual Equation>

[0029] Assuming that a physical quantity desired to obtain is U and thata linear differential operator is A, a nonlinear differential operatoris B and an inhomogeneous term (source term) is f in a nonlinear partialdifferential equation to be satisfied by the physical quantity U, anequation to be solved is expressed as follows:

A−U+B(U)=f  (1)

[0030] This equation can be generally solved by linearly approximating anonlinear term B(U) and by using the successive approximation such asthe SOR method or the ADI method.

[0031] In this embodiment, the following solution method is employed:

[0032] The solution U of equation (1) is represented by using anapproximate solution U^(m) and a perturbation quantity φ (namely, adifference from a true solution) as follows:

U=U ^(m)+φ  (2)

[0033] Also, a nonlinear residual r^(m) of the approximate solutionU^(m) is defined as follows:

r ^(m) =f−A·U ^(m) −B(U ^(m))  (3)

[0034] In the solution method of this embodiment, the perturbationquantity φ is obtained and the approximate solution U^(m) is correctedso as to reduce a norm of a nonlinear residual r^(m+1) employed in asubsequent step, thereby obtaining the true solution U of equation (1).

[0035] When the nonlinear term B(U) is represented by using a nonlinearcoefficient and a differential operator as follows:${B(U)} = {\sum\limits_{n}{{\overset{n}{N}(U)}{\overset{n}{D} \cdot U}}}$

${\overset{n}{N}(U)}\quad \cdots$

[0036] . . . nonlinear coefficient, ${\overset{n}{D}(U)}\quad \cdots$

[0037] . . . differential operator N(U^(m)+φ) is expanded by the Taylorexpansion or the mean value theorem, so as to be expressed, with termshigher than φ² ignored, as follows: $\begin{matrix}{{{\left. {{B\left( U^{m} \right)} + \varphi} \right) \cong {{B\left( U^{m} \right)} + {L \cdot \varphi} + {\varphi {\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot \varphi}}}}}},\quad \ldots}{{{L \cdot \varphi} \equiv {\sum\limits_{n}{\left\lbrack {{{\overset{n}{N}\left( U^{m} \right)}\overset{n}{D}} + {{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot U^{m}}}} \right\rbrack \varphi}}},{{\overset{n}{N^{\prime}}\left( U^{m} \right)} \equiv \left\lbrack \frac{\partial{\overset{n}{N}(U)}}{\partial U} \right\rbrack_{U = U^{m}}}}} & (4)\end{matrix}$

[0038] On the basis of equations (1) through (4), the following holds:

A·(U ^(m)+φ)+B(U ^(m)+φ)=f$\left. \Rightarrow{{A \cdot \left( {U^{m} + \varphi} \right)} + \left( {{B\left( U^{m} \right)} + {L \cdot \varphi} + {\varphi {\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot \varphi}}}}} \right)} \right. = {\left. f\Rightarrow{\left\lbrack {A + L + {\varphi {\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}\overset{n}{D}}}}} \right\rbrack \cdot \varphi} \right. = {{f - {A \cdot U^{m}} - {B\left( U^{m} \right)}} = r^{m}}}$

[0039] Accordingly, in order to obtain the perturbation quantity A, thefollowing equation is solved: $\begin{matrix}{{\left\lbrack {A + L + {\varphi {\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}\overset{n}{D}}}}} \right\rbrack \cdot \varphi} = r^{m}} & (5)\end{matrix}$

[0040] <Algorithm for Nonlinear Residual Cutting Method>

[0041] In the algorithm of this invention, a convergence solution ofequation (5) is not to be obtained but, for example, from the followingterm on the left-hand side of equation (5):$\varphi^{m}{\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot \varphi^{m}}}}$

[0042] φ^(m) is transposed in the form of φ^(m) _(old) to the right-handside, so as to obtain the following linearized equation:${\left\lbrack {A + L} \right\rbrack \cdot \varphi} = {r^{m} - {\varphi_{old}^{m}{\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot \varphi_{old}^{m}}}}}}$

[0043] This equation is solved by the ADI method or the like, so as toobtain a predicted approximate value φ^(m) through a minimum unit ofrepeated calculations with the steepest convergence gradient, and theobtained predicted approximate value φ^(m) is supplied to anoptimization control routine. A series of operations are repeatedlyexecuted until the execution result of the optimization control routinesatisfies a predetermined condition.

[0044] <Optimization Control Routine>

[0045] Assuming that the number of repeating times is m, a synthesizeperturbation quantity φ^(m) for minimizing the norm of the nonlinearresidual and a new approximate solution U^(m+1) are defined as follows:$\begin{matrix}{\varphi^{m} = {{\alpha_{1}\psi^{m}} + {\sum\limits_{l = 2}^{\quad {L\quad \max}}\quad {\alpha_{l}\varphi^{m - l + 1}}}}} & (6)\end{matrix}$

 U^(m+1) =U ^(m)+φ^(m)  (7)

[0046] In equation (6), α₁(1=1, 2, 3, . . . , L_(max)) is a residualminimization coefficient, which is a constant obtained by a calculationmethod described later. The following minimization of the nonlinearresidual is carried out so as to make the norm of the nonlinear residualsmaller:

[0047] When the new approximate solution U^(m+1) is represented byequation (7), the nonlinear residual r^(m+1) of this approximatesolution U^(m+1) is represented, on the basis of equation (3), asfollows:

r ^(m+1) =f−A·U ^(m+1) −B(U ^(m+1))  (8)

[0048] Equation (8) can be expressed by using equations (3) and (4) asfollows: $\begin{matrix}{{r^{m + 1} \cong {f - {A \cdot \left( {U^{m} + \varphi^{m}} \right)} - {B\left( U^{m} \right)} - {L \cdot \varphi^{m}} - {\varphi^{m}{\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot \varphi^{m}}}}}}}\quad = {r^{m} - {\left\lbrack {A + L + {\varphi^{m}{\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}\overset{n}{D}}}}} \right\rbrack \cdot \varphi^{m}}}} & (9)\end{matrix}$

[0049] At this point, when

Ã=A+L  (10)

[0050] equation (9) can be expressed as follows: $\begin{matrix}{r^{m + 1} = {r^{m} - {\varphi^{m}{\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot \varphi^{m}}}}} - {\overset{\sim}{A} \cdot \varphi^{m}}}} & (11)\end{matrix}$

[0051] When, for example, the L² norm is used as the norm for evaluatingthe amplitude of the nonlinear residual, a norm ∥r^(m+1)∥ of thenonlinear residual r^(m+1) is, in the three-dimensional case, given bythe following equation as a square root of a sum of all the inner pointsof (r^(m+1))²: $\begin{matrix}{{r^{m + 1}} = \sqrt{\sum\limits_{ijk}\left( {r^{m} - {\varphi^{m}{\sum\limits_{n}{\overset{n}{N^{\prime}}\left( U^{m} \right){\overset{n}{D} \cdot \varphi^{m}}}}} - {\overset{\sim}{A} \cdot \varphi^{m}}} \right)^{2}}} & (12)\end{matrix}$

[0052] wherein i, j, and k represent a three-dimensional coordinate.

[0053] The residual minimization coefficient α₁(1=1, 2, 3, . . . ,L_(max)) is determined so as to minimize the norm ∥r^(m+1) ∥ of thenonlinear residual r^(m+1) by the least squares method. Specifically,the following equation is to be satisfied: $\begin{matrix}{\frac{\partial{r^{m + 1}}^{2}}{\partial\alpha_{l}} = {0\quad \left( {{l = 1},2,3,\cdots \quad,L_{\max}} \right)}} & (13)\end{matrix}$

[0054] Since the following term:$\varphi^{m}{\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot \varphi^{m}}}}$

[0055] is present on the right-hand side of equation (12), equation (13)is nonlinear simultaneous equations with respect to α₁, and when thisequation is numerically solved through repeated calculations throughlinear approximation, α₁ can be defined. When α₁ is defined, φ^(m) isfound on the basis of equation (6), and U^(m+1) is found on the basis ofequation (7). The method for calculating the residual minimizationcoefficient α₁ will be described later.

[0056] While incrementing the number m, similar routines are repeated,so as to attain convergence of the solution U for making zero orminimizing the norm of the nonlinear residual of equation (8).

[0057] Now, the control method of this embodiment will be described withreference to the flowchart of FIG. 2.

[0058] First, in step SA1, an initial value U⁰ of the physical quantityU to be obtained is set. In step SA2, the number m of repeating times isset to an initial value of 0, and an initial value r⁰ of the nonlinearresidual r is set to (f−A·U⁰−B(U⁰)).

[0059] Thereafter, procedures of steps SA3 through SA14 are repeatedlyexecuted with the number m of repeating times incremented until theapproximate solution U^(m) is converged.

[0060] In steps SA3 through SA7, the approximate solution (predictedapproximate value) φ^(m) of the following equation is obtained by usingan internal solver that executes the successive approximation such asthe ADI method or the SOR method by linearizing this equation:${\left\lbrack {\overset{\sim}{A} + {\varphi {\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}\overset{n}{D}}}}} \right\rbrack \cdot \varphi} = r^{m}$

[0061] In this embodiment, the upper limit N_(max) of the number ofrepeating times of the successive calculations is set (steps SA5 andSA6).

[0062] When the predicted approximate value φ⁰ is obtained through theprocedures of steps SA3 through SA7, the optimization routine is carriedout by using this predicted approximate value φ⁰ so as to obtain acorrected value φ⁰ for minimizing a nonlinear residual r¹ of asubsequent step. By using this corrected value φ⁰, a primary approximatevalue U¹ of the physical quantity and a primary nonlinear residual r¹can be obtained as follows on the basis of equations (2) and (3):

U ¹ =U ⁰+φ⁰  (14)

r ¹ =f−A·U ¹ −B(U ¹)  (15)

[0063] In step SA13, the number m of repeating times is incremented, andthe flow returns to step SA3, so as to repeat similar calculations. Byrepeating these procedures, (U², r²), (U³, r³), . . . , (U^(m), r^(m)) .. . can be successively obtained, and the norm ∥r^(m)∥ of the nonlinearresidual is reduced.

[0064] When the number of repeating times is m, an equation used forobtaining the predicted approximate value φ^(m) in step SA4 is asfollows:${\left\lbrack {\overset{\sim}{A} + {\varphi {\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}\overset{n}{D}}}}} \right\rbrack \cdot \varphi} = r^{m}$

[0065] Therefore, on the basis of the corrected approximate value φ^(m)having been optimized through the optimization control routine performedin steps SA8 through SA11, the following equations are obtained:

U^(m+1) =U ^(m)+φ^(m)

r^(m+1) f−A·U ^(m+1) −B(U ^(m+1))

[0066] <Method for Obtaining Residual Minimization Coefficient α₁>

[0067] At this point, for simplifying the expression, δ is introduced asfollows:

δ^(m)=φ^(m)

δ^(m−l+1)=φ^(m−l+1)(wherein l=2, 3, . . . , L_(max))

[0068] In this case, the corrected approximate value φ^(m) is expressedas follows: $\begin{matrix}{\varphi^{m} = {\sum\limits_{l = 1}^{L\quad \max}\quad {\alpha_{l}\delta^{m - l + 1}}}} & (17)\end{matrix}$

[0069] Therefore, the following equation holds: $\begin{matrix}\begin{matrix}{r^{m + 1} = {r^{m} - {A \cdot \varphi^{m}} - {B\left( \varphi^{m} \right)}}} \\{= {r^{m} - {\overset{\sim}{A} \cdot \varphi^{m}} - {\varphi^{m}{\sum\limits_{n}\quad {{\overset{n}{N}'}\left( U^{m} \right){\overset{n}{D} \cdot \varphi^{m}}}}}}}\end{matrix} & (18)\end{matrix}$

[0070] When this equation is substituted in equation (17), the followingis obtained: $\begin{matrix}\begin{matrix}{r^{m + 1} = {r^{m} - {\sum\limits_{l = 1}^{L_{\max}}\quad {\alpha_{l}{\overset{\sim}{A} \cdot \delta^{m - l + 1}}}} -}} \\{{\sum\limits_{l = 1}^{L_{\max}}\quad {\sum\limits_{{l'} = 1}^{L_{\max}}\left( {\alpha_{l}\alpha_{l'}\delta^{m - {l'} + 1}{\sum\limits_{n}\quad {{\overset{n}{N}'}\left( U^{m} \right){\overset{n}{D} \cdot \delta^{m - {l'} + 1}}}}} \right)}}}\end{matrix} & (19)\end{matrix}$

[0071] The residual minimization coefficient α₁ is obtained by using theleast squares method under a condition that the L² norm of the nonlinearresidual of the (m+1)th order is to be minimized. $\begin{matrix}\begin{matrix}{{r^{m + 1}}^{2} = {\sum\limits_{ijk}\quad \left( {r^{m} - {\sum\limits_{l = 1}^{L_{\max}}\quad {\alpha_{l}{\overset{\sim}{A} \cdot \delta^{m - l + 1}}}} -} \right.}} \\\left. {\sum\limits_{l = 1}^{L_{\max}}\quad {\sum\limits_{{l'} = 1}^{L_{\max}}\left( {\alpha_{l}\alpha_{l'}\delta^{m - {l'} + 1}{\sum\limits_{n}\quad {{\overset{n}{N}'}\left( U^{m} \right){\overset{n}{D} \cdot \delta^{m - {l'} + 1}}}}} \right)}} \right)^{2} \\{\square}\end{matrix} & (20)\end{matrix}$

[0072] When this equation is differentiated with respect to α₁, thefollowing is obtained: $\begin{matrix}{\frac{\partial{r^{m + 1}}^{2}}{\partial\alpha_{l}} = {0\quad \left( {{l = 1},2,3,\ldots \quad,L_{\max}} \right)}} & (21)\end{matrix}$

[0073] When this equation (21) is solved with respect to α₁, theresidual minimization coefficient α₁ can be obtained. This equation isgenerally nonlinear simultaneous equations with respect to α₁ and issolved through linearization and repeated calculations. For example,when φ^(m) of the third term on the right-hand side of equation (18) istransposed in the form of φ^(m) _(old) to the right-hand side to belinearized and the following is defined: $\begin{matrix}{{\overset{\sim}{r}}^{m} = {r^{m} - {\varphi_{old}^{m}{\sum\limits_{n}\quad {{\overset{n}{N}'}\left( U^{m} \right){\overset{n}{D} \cdot \varphi_{old}^{m}}}}}}} & (22)\end{matrix}$

[0074] the following equation holds: $\begin{matrix}{{r^{m + 1} = {{\overset{\sim}{r}}^{m} - {\sum\limits_{l = 1}^{L_{\max}}\quad {\alpha_{l}{\overset{\sim}{A} \cdot \delta^{m - l + 1}}}}}}{{r^{m + 1}}^{2} = {\sum\limits_{ijk}\quad \left( {{\overset{\sim}{r}}^{m} - {\sum\limits_{l = 1}^{L_{\max}}\quad {\alpha_{l}{\overset{\sim}{A} \cdot \delta^{m - l + 1}}}}} \right)^{2}}}} & (23)\end{matrix}$

[0075] When this equation is substituted in equation (21) and theresultant is differentiated with respect to α₁, the following isobtained: $\begin{matrix}{{\frac{\partial{r^{m + 1}}^{2}}{\partial\alpha_{l^{\prime}}} = {{- 2}{\sum\limits_{ijk}{\left\lbrack {{\overset{\sim}{r}}^{m} - {\overset{\sim}{A} \cdot {\sum\limits_{l = 1}^{L_{\max}}{\alpha_{l}\delta^{m - l + 1}}}}} \right\rbrack {\overset{\sim}{A} \cdot \delta^{m - l^{\prime} + 1}}}}}},{{\left( {{l^{\prime} = 1},2,3,\ldots \quad,L_{\max}} \right)\therefore{\sum\limits_{ijk}{\sum\limits_{l = 1}^{L_{\max}}{{\alpha_{l}\left( {\overset{\sim}{A} \cdot \varphi^{m - l + 1}} \right)}\left( {\overset{\sim}{A} \cdot \varphi^{m - l^{\prime} + 1}} \right)}}}} = {\sum\limits_{ijk}{{\overset{\sim}{r}}^{m}\left( {\overset{\sim}{A} \cdot \varphi^{m - l^{\prime} + 1}} \right)}}},\left( {{l^{\prime} = 1},2,3,\ldots \quad,L_{\max}} \right)} & (24)\end{matrix}$

[0076] When the L_(max)-element simultaneous linear equations given asequation (24) are solved, α₁ is temporarily obtained, and φ^(m) obtainedby substituting the temporal α₁ in equation (17) is substituted inequation (22) in the form of φ^(m) _(old) again, and the resultantsimultaneous linear equations given as equation (24) are solved again.Such repeated calculations are executed several times, so as to obtainthe residual minimization coefficient α₁.

[0077] Embodiment 2

[0078]FIG. 3 is a flowchart for showing procedures in a control methodaccording to Embodiment 2 of the invention.

[0079] <Formulation of Nonlinear Residual Equation>

[0080] Assuming that a physical quantity desired to obtain is U and thata linear differential operator is A, a nonlinear differential operatoris B and an inhomogeneous term (source term) is f in a nonlinear partialdifferential equation to be satisfied by the physical quantity U, anequation to be solved is, as in Embodiment 1, as follows:

A·U+B(U)=f  (1)

[0081] This equation can be generally solved by linearly approximating anonlinear term B(U) and by using the successive approximation such asthe SOR method or the ADI method.

[0082] In this embodiment, the following solution method is employed:

[0083] The solution of equation (1) is represented by using anapproximate solution U^(m) and a perturbation quantity 0 (namely, adifference from a true solution) as follows:

U=U ^(m)+φ  (2)

[0084] Also, a nonlinear residual r^(m) of the approximate solutionU^(m) is defined as follows:

r ^(m) f−A·U ^(m) −B(U ^(m))  (3)

[0085] In the solution method of this embodiment, the perturbationquantity φ is obtained and the approximate solution U^(m) is correctedso as to reduce a nonlinear residual r^(m+1) employed in a subsequentstep, thereby obtaining the true solution U of equation (1).

[0086] When the nonlinear term B(U) is represented by using a nonlinearcoefficient and a differential operator as follows:${B(U)} = {\sum\limits_{n}\quad {{\overset{n}{N}(U)}{\overset{n}{D} \cdot U}}}$

${\overset{n}{N}(U)}\quad \cdots$

[0087] . . . nonlinear coefficient, ${\overset{n}{D}(U)}\quad \cdots$

[0088] . . . differential operator N(U^(m)+φ) is expanded by the Taylorexpansion or the mean value theorem, so as to be expressed, with termshigher than φ ignored, as follows:

B(U ^(m)+φ)≅B(U ^(m))+L·φ,  (4′)${{L \cdot \varphi} \equiv {\sum\limits_{n}\quad {\left\lbrack {{{\overset{n}{N}\left( U^{m} \right)}\overset{n}{D}} + {{\overset{n}{N}'}\left( U^{m} \right){\overset{n}{D} \cdot U^{m}}}} \right\rbrack \varphi}}},{{{\overset{n}{N}'}\left( U^{m} \right)} \equiv \left\lbrack \frac{\partial{\overset{n}{N}(U)}}{\partial U} \right\rbrack_{U = U^{m}}}$

[0089] On the basis of equations (1) through (4′), the following holds:

A·(U ^(m)+φ)+B(U ^(m)+φ)=f

→A·(U ^(m)+φ)+(B(U ^(m))+L·φ)=f

→[A+L]·φ=f−A·U ^(m) −B(U ^(m))=r ^(m)

[0090] Accordingly, in order to obtain the perturbation quantity φ, thefollowing equation is solved:

[A+L]·φ=r ^(m)  (5′)

[0091] <Algorithm for Nonlinear Residual Cutting Method>

[0092] In the algorithm of this invention, a convergence solution ofequation (5′) is not to be obtained but an equation obtained bylinearizing the left-hand side of equation (5′) is solved by the ADImethod or the like, so as to obtain a predicted approximate value φ^(m)through a minimum unit of repeated calculations with the steepestconvergence gradient, and the obtained predicted approximate value φ^(m)is supplied to an optimization control routine. A series of operationsare repeatedly executed until the execution result of the optimizationcontrol routine satisfies a predetermined condition.

[0093] <Optimization Control Routine>

[0094] Assuming that the number of repeating times is m, a synthesizeperturbation quantity φ^(m) for minimizing the norm of the nonlinearresidual and a new approximate solution U^(m+1) are defined as follows:$\begin{matrix}{\varphi^{m} = {{\alpha_{1}\psi^{m}} + {\sum\limits_{l = 2}^{L_{\max}}\quad {\alpha_{l}\varphi^{m - l + 1}}}}} & (6)\end{matrix}$

 U ^(m+1) =U ^(m)+φ^(m)  (7)

[0095] In equation (6), α₁(1=1, 2, 3, . . . , L_(max)) is a residualminimization coefficient, which is a constant obtained by a calculationmethod described later. The following minimization of the nonlinearresidual is carried out so as to make the norm of the nonlinear residualsmaller:

[0096] When the new approximate solution U^(m+1) is represented byequation (7), the nonlinear residual r^(m+1) of this approximatesolution U^(m+1) is expressed, on the basis of equation (3), as follows:

r ^(m+1) =f−A·U ^(m+1) −B(U ^(m+1))  (8)

[0097] Equation (8) can be expressed by using equations (3) and (4′) asfollows: $\begin{matrix}{{r^{m + 1} \approx {f - {A \cdot \left( {U^{m} + \varphi^{m}} \right)} - {B\left( U^{m} \right)} - {L \cdot \varphi^{m}}}}\quad = {r^{m} - {\left\lbrack {A + L} \right\rbrack \cdot \varphi^{m}}}} & \left( {9'} \right)\end{matrix}$

[0098] At this point, when

Ã=A+L  (10)

[0099] equation (9′) can be expressed as follows:

r ^(m+1) =r ^(m) −Ã·φ^(m)  (11′)

[0100] When, for example, the L² norm is used as the norm for evaluatingthe amplitude of the nonlinear residual, a norm ∥r^(m+1)∥ of thenonlinear residual r^(m+1) is, in the three-dimensional case, given bythe following equation as a square root of a sum of all the inner pointsof (r^(m+1))²: $\begin{matrix}{{r^{m + 1}} = \sqrt{\sum\limits_{ijk}\left( {r^{m} - {\overset{\sim}{A} \cdot \varphi^{m}}} \right)^{2}}} & \left( {12'} \right)\end{matrix}$

[0101] The residual minimization coefficient α₁(1=1, 2, 3, . . . ,L_(max)) is determined so as to minimize the norm ∥r^(m+1)∥ of thenonlinear residual r_(m+1) by the least squares method. Specifically,the following equation is to be satisfied: $\begin{matrix}{\frac{\partial{r^{m + 1}}^{2}}{\partial\alpha_{l}} = {0\quad \left( {{l = 1},2,3,\cdots \quad,L_{\max}} \right)}} & (13)\end{matrix}$

[0102] Equation (13) is the L_(max)-element simultaneous equations, andwhen this equation is numerically solved, α₁ can be defined. When α₁ isdefined, φ^(m) is found on the basis of equation (6), and U^(m+1) isfound on the basis of equation (7). The method for calculating theresidual minimization coefficient α₁ will be described later.

[0103] While incrementing the number m, similar routines are repeated,so as to attain convergence of the solution U for making zero orminimizing the norm of the nonlinear residual of equation (8).

[0104] Now, the control method of this embodiment will be described withreference to the flowchart of FIG. 3.

[0105] First, in step SB1, an initial value U⁰ of the physical quantityU to be controlled is set. In step SB2, the number m of repeating timesis set to an initial value of 0, and an initial value r⁰ of thenonlinear residual r is set to (f−A·U⁰−B(U⁰)).

[0106] Thereafter, procedures of steps SB3 through SB12 are repeatedlyexecuted with the number m of repeating times incremented until theapproximate solution U^(m) is converged.

[0107] In steps SB3 through SB6, the approximate solution (predictedapproximate value) φ^(m) of the following equation is obtained by usingan internal solver that executes the existing successive approximationsuch as the ADI method or the SOR method:

Ã·φ=r ^(m)

[0108] However, when this equation is to be solved by the existingsuccessive approximation, a very large number of repeated times isrequired before the convergence for practical use, and hence thecalculation takes massive time.

[0109] Therefore, in this embodiment, the upper limit N_(max) of thenumber of repeating times of the successive calculations is set (stepsSB5 and SB6).

[0110] When the predicted approximate value φ⁰ is obtained through theprocedures of steps SB3 through SB6, the optimization routine is carriedout, in steps SB8 and SB9, by using this predicted approximate value φ⁰so as to obtain a corrected value φ⁰ for minimizing a subsequentnonlinear residual r¹. By using this corrected value φ⁰, a primaryapproximate value U¹ of the physical quantity and a primary nonlinearresidual r¹ can be obtained as follows on the basis of equations (2) and(3):

U ¹ =U ⁰+φ⁰  (14)

r ¹ =f−A·U ¹ −B(U ¹)  (15)

[0111] In step SB13, the number m of repeating times is incremented, andthe flow returns to step SB3, so as to repeat similar calculations. Byrepeating these procedures, (U², r²), (U³, r³), . . . , (U^(m), r^(m)) .. . can be successively obtained, and the norm ∥r^(m)∥ of the nonlinearresidual is reduced.

[0112] When the number of repeating times is m, an equation used forfinding the predicted approximate value φ^(m) in step SB4 is as follows:

Ã·φ=r ^(m)

[0113] Therefore, on the basis of the corrected approximate value φ^(m)having been optimized through the optimization control routine performedin steps SB8 and SB9, the following equations are obtained in step SB12:

U ^(m+1) =U ^(m)+φ^(m)

r ^(m+1) =f−A·U ^(m+1) −B(U ^(m+1))

[0114] <Method for Obtaining Residual Minimization Coefficient α₁>

[0115] At this point, for simplifying the expression, δ is introduced asfollows:

δ^(m)=φ^(m)

δ^(m−l+1)=φ^(m−l+1) (wherein l=2, 3, . . . , L_(max))

[0116] In this case, the corrected approximate value φ^(m) is expressedas follows: $\begin{matrix}{\varphi^{m} = {\sum\limits_{l = 1}^{L\quad \max}\quad {\alpha_{l}\delta^{m - l + 1}}}} & (17)\end{matrix}$

[0117] Therefore, the following equation holds: $\begin{matrix}{{r^{m + 1} = {{r^{m} - {A \cdot \varphi^{m}} - {B\left( \varphi^{m} \right)}}\quad = {r^{m} - {\overset{\sim}{A} \cdot \varphi^{m}}}}}\quad} & \left( {18'} \right)\end{matrix}$

[0118] When this equation is substituted in equation (17), the followingis obtained: $\begin{matrix}{r^{m + 1} = {r^{m} - {\sum\limits_{l = 1}^{L\quad \max}\quad {\alpha_{l}{\overset{\sim}{A} \cdot \delta^{m - l + 1}}}}}} & \left( {19'} \right)\end{matrix}$

[0119] The residual minimization coefficient α₁ is obtained by using theleast squares method under a condition that the L² norm of the nonlinearresidual of the (m+1)th order is to be minimized. $\begin{matrix}{{r^{m + 1}}^{2} = {\sum\limits_{ijk}\left( {r^{m} - {\sum\limits_{l = 1}^{L\quad \max}\quad {\alpha_{l}{\overset{\sim}{A} \cdot \delta^{m - l + 1}}}}} \right)^{2}}} & \left( {20'} \right)\end{matrix}$

[0120] When this equation is differentiated with respect to α₁, thefollowing is obtained: $\begin{matrix}{\frac{\partial{r^{m + 1}}^{2}}{\partial\alpha_{l}} = {0\quad \left( {{l = 1},2,3,\cdots \quad,L_{\max}} \right)}} & (21)\end{matrix}$

[0121] As a result of differentiation with respect to α₁, the followingis obtained: $\begin{matrix}{{\frac{\partial{r^{m + 1}}^{2}}{\partial\alpha_{r}} = {{- 2}{\sum\limits_{ijk}{\left\lbrack {r^{m} - {\overset{\sim}{A} \cdot {\sum\limits_{l = 1}^{L\quad \max}\quad {\alpha_{l}\delta^{m - l + 1}}}}} \right\rbrack {\overset{\sim}{A} \cdot \delta^{m - l + 1}}}}}},{{\left( {{l^{\prime} = 1},2,3,\cdots \quad,L_{\max}} \right)\therefore{\sum\limits_{ijk}{\sum\limits_{l = 1}^{L\quad \max}\quad {{\alpha_{l}\left( {\overset{\sim}{A} \cdot \varphi^{m - l + 1}} \right)}\left( {\overset{\sim}{A} \cdot \varphi^{m - l^{\prime} + 1}} \right)}}}} = {\sum\limits_{ijk}{r^{m}\left( {\overset{\sim}{A} \cdot \varphi^{m - l^{\prime} + 1}} \right)}}},\left( {{l^{\prime} = 1},2,3,\cdots \quad,L_{\max}} \right)} & (25)\end{matrix}$

[0122] When the L_(max)-element simultaneous linear equations given asequation (25) are solved, the residual minimization coefficient α₁ canbe obtained.

[0123] In the case where the algorithm according to Embodiment 1 or 2 isapplied to the drive control for a vehicle shown in FIG. 1, a predictedcorrection calculation unit 14 first obtains, through repeatedcalculations, predicted correction φ^(m) on the basis of the currentdrive information. Next, an optimal addition correction calculation unit15 calculates an addition correction φ^(m) of the mth-order. Acorrection convergence judgment unit 16 ultimately determines whether ornot the solutions obtained by the predicted correction calculation unit14 and the optimal addition correction calculation unit 15 have beenconverged, so as to determine the optimal control for returning thevehicle to the predetermined orbit. The determined control informationis sent to a direction control unit 17 for correcting the drive route ofthe vehicle.

[0124] Although the drive control for a vehicle is herein exemplified todescribe the invention, it goes without saying that the invention isapplicable to control employed in any of various fields.

[0125] Also, the control method can be easily realized by an apparatusequipped with a computer for executing a program for realizing thecontrol method. Furthermore, when the program for realizing the controlmethod is recorded in a computer-readable recording medium, the controlmethod can be realized by allowing a computer to execute the programrecorded in the recording medium.

[0126] In approximately expressing the nonlinear term, B(U^(m)+φ), termsof higher orders than the second-order are ignored in the expansionequation of the nonlinear coefficient in Embodiment 1 and terms ofhigher orders than the first-order are ignored in Embodiment 2, but theorder in which the expansion is ended is not limited to them. However,in consideration of programmability and the like, it seems that theexpansion is practically ended in the first-order or the second-order.

[0127] <Application to Numerical Calculation>

[0128] The algorithm according to this invention is also applicable tonumerical calculation.

[0129]FIGS. 4A and 4B show, as an example of the application tonumerical calculation, application to nonlinear simultaneous equationsof a two-dimensional Thompson grid generation method. Also, FIGS. 5A and5B are graphs for showing comparison in the convergence rate in theexemplified application shown in FIGS. 4A and 4B, and specifically, FIG.5A is a graph obtained by utilizing a conventional linear residualcutting method and FIG. 5B is a graph obtained by using the nonlinearresidual cutting method according to Embodiment 1 of the invention. Inthese graphs, the ordinate indicates the norm of the nonlinear residual(expressed in logarithm) and the abscissa indicates the number ofrepeating times. In each graph, the norm of the nonlinear residual of anequation for finding the x-coordinate of a grid point is plotted byusing a thin line, and the norm of the nonlinear residual of an equationfor finding the y-coordinate of a grid point is plotted by using a thickline. It is understood from FIGS. 5A and 5B that the convergence rate isremarkably improved by the invention.

[0130] Similarly, FIGS. 6A and 6B show, as an example of the applicationto numerical calculation, application to the nonlinear simultaneousequations of the two-dimensional Thompson grid generation method. Also,FIGS. 7A and 7B are graphs for showing comparison in the convergencerate in the exemplified application shown in FIGS. 6A and 6B, andspecifically, FIG. 7A is a graph obtained by utilizing the conventionallinear residual cutting method and FIG. 7B is a graph obtained by usingthe nonlinear residual cutting method according to Embodiment 2 of theinvention. In these graphs, the ordinate indicates the norm of thenonlinear residual (expressed in logarithm) and the abscissa indicatesthe number of repeating times. In each graph, the norm of the nonlinearresidual of an equation for finding the x-coordinate of a grid point isplotted by using a thin line, and the norm of the nonlinear residual ofan equation for finding the y-coordinate of a grid point is plotted byusing a thick line. It is understood from FIGS. 7A and 7B that even whena solution is not converged by the conventional method, it can beconverged by the method of this invention.

[0131] Such application to numerical calculation can be easily realizedby an apparatus equipped with a computer for executing a program forrealizing the calculation. Furthermore, when the program for realizingthe calculation is recorded in a computer-readable recording medium, thecalculation can be realized by allowing a computer to execute theprogram recorded in the recording medium.

[0132] As described so far, since the present invention employs thesuccessive approximation algorithm in which a nonlinear residual isdirectly dealt with so as to be reduced, definite and rapid control fora nonlinear operation with strong nonlinearity can be easily performed.

What is claimed is:
 1. A method for controlling a physical quantity U,comprising: a step of solving, by successive approximation, A·U+B(U)=f,wherein A is a linear differential operator, B is a nonlineardifferential operator, and f is an inhomogeneous term (source term) in anonlinear partial differential equation to be satisfied by said physicalquantity U, wherein (f−A·U^(m)−B(U^(m))) is given as a nonlinearresidual r^(m) of an approximate solution U^(m), wherein m is the numberof repeating times, and said approximate solution U^(m) is corrected soas to reduce a nonlinear residual r^(m+1) employed in a subsequent step.2. The method for controlling a physical quantity U of claim 1, whereinthe nonlinear term B(U) is expressed as follows:${B(U)} = {\sum\limits_{n}{{\overset{n}{N}(U)}{\overset{n}{D} \cdot U}}}$

${\overset{n}{N}(U)}\quad \cdots$

. . . nonlinear coefficient, ${\overset{n}{D}(U)}\quad \cdots$

. . . differential operator N(U^(m)+φ) is expanded with respect to saidapproximate solution U^(m) and a perturbation quantity φ by the Taylorexpansion or the mean value theorem, to give the following with termshigher than φ² ignored: $\begin{matrix}{{{B\left( {U^{m} + \varphi} \right)} \cong {{B\left( U^{m} \right)} + {L \cdot \varphi} + {\varphi \quad {\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot \varphi}}}}}},{{L \cdot \varphi} \equiv {\sum\limits_{n}{\left\lbrack {{{\overset{n}{N}\left( U^{m} \right)}\overset{n}{D}} + {{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot U^{m}}}} \right\rbrack \varphi}}},{{\overset{n}{N^{\prime}}\left( U^{m} \right)} \equiv \left\lbrack \frac{\partial{\overset{n}{N}(U)}}{\partial U} \right\rbrack_{U = U^{m}}}} & (4)\end{matrix}$

wherein a first process and a second process are repeatedly executeduntil said approximate solution U^(m) is converged, wherein said firstprocess includes the steps of: setting U⁰ as an initial value of saidphysical quantity U; setting 0 as an initial value of said number m ofrepeating times and (f−A·U⁰−B(U⁰)) as an initial value r⁰ of a nonlinearresidual r; and obtaining a predicted approximate value φ^(m) of thefollowing equation through repeated calculations while incrementing saidnumber m of repeating times: $\begin{matrix}{{\left\lbrack {A + L + {\varphi \quad {\sum\limits_{n}{{\overset{n}{N^{\prime}}\left( U^{m} \right)}\overset{n}{D}}}}} \right\rbrack \cdot \varphi} = r^{m}} & (5)\end{matrix}$

and wherein said second process includes the steps of: obtaining acorrected approximate value φ^(m) for minimizing a norm of a nonlinearresidual r^(m+1) on the basis of a predicted approximate value φ^(m) andcorrected approximate values φ^(m−1), . . . , and φ^(m−Lmax+1), whereinLmax is an integer of 2 or more; giving (U^(m)+φ^(m)) as an approximatesolution U^(m+1); and giving (f−A·U^(m+1)−B(U^(m+1))) as said nonlinearresidual r^(m+1).
 3. The method for controlling a physical quantity U ofclaim 1, wherein the nonlinear term B(U) is expressed as follows:${B(U)} = {\sum\limits_{n}{{\overset{n}{N}(U)}{\overset{n}{D} \cdot U}}}$

${\overset{n}{N}(U)}\quad \cdots$

. . . nonlinear coefficient, ${\overset{n}{D}(U)}\quad \cdots$

. . . differential operator N(U^(m)+φ) is expanded with respect to saidapproximate solution U^(m) and a perturbation quantity φ by the Taylorexpansion or the mean value theorem, to give the following with termshigher than φ ignored: B(U ^(m)+φ)≅B(U ^(m))+L·φ,  (4′)${{L \cdot \varphi} \equiv {\sum\limits_{n}{\left\lbrack {{{\overset{n}{N}\left( U^{m} \right)}\overset{n}{D}} + {{\overset{n}{N^{\prime}}\left( U^{m} \right)}{\overset{n}{D} \cdot U^{m}}}} \right\rbrack \varphi}}},{{\overset{n}{N^{\prime}}\left( U^{m} \right)} \equiv \left\lbrack \frac{\partial{\overset{n}{N}(U)}}{\partial U} \right\rbrack_{U = U^{m}}}$

wherein a first process and a second process are repeatedly executeduntil said approximate solution U^(m) is converged, wherein said firstprocess includes the steps of: setting U⁰ as an initial value of saidphysical quantity U; setting 0 as an initial value of said number m ofrepeating times and (f−A·U⁰−B(U⁰)) as an initial value r⁰ of a nonlinearresidual r; and obtaining a predicted approximate value φ^(m) of[A+L]φ=r^(m) through repeated calculations while incrementing saidnumber m of repeating times, and wherein said second process includesthe steps of: obtaining a corrected approximate value φ^(m) forminimizing a norm of a nonlinear residual r^(m+1) on the basis of apredicted approximate value φ^(m) and corrected approximate valuesφ^(m−1), . . . , and φ^(m−Lmax+1), wherein Lmax is an integer of 2 ormore; giving (U^(m)+φ) as an approximate solution U^(m+1); and giving(f−A·U^(m+1)−B(U^(m+1))) as said nonlinear residual r^(m+1).
 4. Acontroller for controlling a physical quantity U, comprising: means forsolving, by successive approximation, A·U+B(U)=f, wherein A is a lineardifferential operator, B is a nonlinear differential operator, and f isan inhomogeneous term (source term) in a nonlinear partial differentialequation to be satisfied by said physical quantity U, wherein(f−A·U^(m)−B(U^(m))) is given as a nonlinear residual r^(m) of anapproximate solution U^(m), wherein m is the number of repeating times,and said approximate solution Um is corrected so as to reduce anonlinear residual r^(m+1) employed in a subsequent step.
 5. A recordingmedium in which a control program for allowing a computer to control aphysical quantity U is recorded, said program allowing said computer toexecute: a processing for solving, by successive approximation,A·U+B(U)=f, wherein A is a linear differential operator, B is anonlinear differential operator, and f is an inhomogeneous term (sourceterm) in a nonlinear partial differential equation to be satisfied bysaid physical quantity U, a processing for giving (f−A·U^(m)−B(U^(m)))as a nonlinear residual r^(m) of an approximate solution U^(m), whereinm is the number of repeating times, and a processing for correcting saidapproximate solution U^(m) so as to reduce a nonlinear residual r^(m+1)employed in a subsequent step.
 6. A method for numerically calculating aphysical quantity U, comprising: a step of solving, by successiveapproximation, A·U+B(U)=f, wherein A is a linear differential operator,B is a nonlinear differential operator, and f is an inhomogeneous term(source term) in a nonlinear partial differential equation to besatisfied by said physical quantity U, wherein (f−A·U^(m)−B(U^(m))) isgiven as a nonlinear residual rr of an approximate solution U^(m),wherein m is the number of repeating times, and said approximatesolution U^(m) is corrected so as to reduce a nonlinear residual r^(m+1)employed in a subsequent step.
 7. A numerical calculator for calculatinga physical quantity U, comprising: means for solving, by successiveapproximation, A·U+B(U)=f, wherein A is a linear differential operator,B is a nonlinear differential operator, and f is an inhomogeneous term(source term) in a nonlinear partial differential equation to besatisfied by said physical quantity U, wherein (f−A·U^(m)−B(U^(m))) isgiven as a nonlinear residual r^(m) of an approximate solution U^(m),wherein m is the number of repeating times, and said approximatesolution U^(m) is corrected so as to reduce a nonlinear residual r^(m+1)employed in a subsequent step.
 8. A recording medium in which anumerical calculation program for allowing a computer to numericallycalculate a physical quantity U, said numerical calculation programallowing said computer to execute: a processing for solving, bysuccessive approximation, A·U+B(U)=f, wherein A is a linear differentialoperator, B is a nonlinear differential operator, and f is aninhomogeneous term (source term) in a nonlinear partial differentialequation to be satisfied by said physical quantity U, a processing forgiving (f−A·U^(m)−B(U^(m))) as a nonlinear residual rr of an approximatesolution U^(m), wherein m is the number of repeating times, and aprocessing for correcting said approximate solution U^(m) so as toreduce a nonlinear residual r^(m+1) employed in a subsequent step.